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ABSTRACT 

The discovery by Marcy et al. (2001) of two planets in 2:1 orbital resonance about the star 
GJ876 has been supplemented by a dynamical fit to the data by Laughlin & Chambers (2001) 
which places the planets in coplanar orbits deep in three resonances at the 2:1 mean-motion 
commensurability. The selection of this almost singular state by the dynamical fit means that the 
resonances are almost certainly real, and with the small amplitudes of libration of the resonance 
variables, indefinitely stable. Several unusual properties of the 2:1 resonances are revealed by 
the GJ 876 system. The libration of both lowest order mean-motion resonance variables and the 
secular resonance variable, 6\ = Ai — 2A2 +roi, 62 = Ai — 2A2 + 1172, and $3 = vj\ — w-i, about 0° 
(where A1.2 are the mean longitudes of the inner and outer planet and 071,2 are the longitudes of 
periapse) differs from the familiar geometry of the Io-Europa pair, where #2 and #3 librate about 
180°. By considering the condition that w\ = W2 for stable simultaneous librations of 9\ and 62, 
we show that the GJ876 geometry results because of the large orbital eccentricities e,, whereas 
the very small eccentricities in the Io-Europa system lead to the latter's geometry. Surprisingly, 
the GJ 876 configuration, with 6\, 62, and #3 all librating, remains stable for e\ up to 0.86 and for 
amplitude of libration of 9\ approaching 45° with the current eccentricities — further supporting 
the indefinite stability of the existing system. 

Any process that drives originally widely separated orbits toward each other could result 
in capture into the observed resonances at the 2:1 commensurability. We find that forced in- 
ward migration of the outer planet of the GJ 876 system results in certain capture into the 
observed resonances if initially e± < 0.06 and e2 < 0.03 and the migration rate ja^/aal ^ 
3 x 10~ 2 (a2/ AU) -3 / 2 yr _1 . Larger eccentricities lead to likely capture into higher order res- 
onances before the 2:1 commensurability is reached. The planets are sufficiently massive to open 
gaps in the nebular disk surrounding the young GJ 876 and to clear the disk material between 
them, and the resulting planet-nebular interaction typically forces the outer planet to migrate 
inward on the disk viscous time scale, whose inverse is about three orders of magnitude less 
than the above upper bound on \a,2/a2\ for certain capture. If there is no eccentricity damping, 
eccentricity growth is rapid with continued migration within the resonance, with exceeding the 
observed values after a further reduction in the semi-major axes of only 7%. With eccentricity 
damping e^/e^ = — K\a,i/ai\, the eccentricities reach equilibrium values that remain constant for 
arbitrarily long migration within the resonances. The equilibrium eccentricities are close to the 
observed eccentricities for K 100 if there is migration and damping of the outer planet only, 
but for K ss 10 if there is also migration and damping of the inner planet. This result is indepen- 
dent of the magnitude or functional form of the migration rate hi as long as e^/e^ = — K\hi/ai\. 
Although existing analytic estimates of the effects of planet-nebula interaction are consistent with 
this form of eccentricity damping for certain disk parameter values, it is as yet unclear that such 
interaction can produce the large value of K required to obtain the observed eccentricities. The 
alternative eccentricity damping by tidal dissipation within the star or the planets is completely 
negligible, so the observed dynamical properties of the GJ 876 system may require an unlikely 
fine tuning of the time of resonance capture to be near the end of the nebula lifetime. 



1. INTRODUCTION 

Marcy et al. (2001) have discovered two planets 
about the nearby M dwarf star GJ876. A prelim- 
inary fit of the stellar radial velocity (RV) vari- 
ations due to two unperturbed Kepler orbits im- 
plies that the orbital periods of the two planets 
are nearly in the ratio 2:1. This resonance is an 
analog to the orbital resonances among the satel- 
lites of Jupiter and Saturn (e.g., Peale 1999), but 
it is the first to be discovered among extrasolar 
planetary systems. Table 1 shows the system pa- 
rameters from the Marcy et al. analysis based on 
data taken at both the Keck and Lick Observato- 
ries and an adopted stellar mass of O.32M . The 
orbital periods are approximately 30 and 60 days. 
The reduced chi-square statistic of xt — 1-88 indi- 
cates that the two-Kepler system is an adequate fit 
to the data. However, the rather large minimum 
masses of the planets of 0.56 and 1.89 Mj (Jupiter 
masses) mean that the two-Kepler fit may not be 
a good determination of the system characteris- 
tics because the large mutual perturbations will 
ensure that the orbits deviate from Kepler orbits. 
In fact, substitution of the Marcy et al. parame- 
ters from Table 1 (with sinz = 1 for both planets, 
where i is the inclination of the orbital plane to 
the plane of the sky) into a calculation of the per- 
turbed orbital motions beginning at the specified 
initial epoch leads to large variations in the or- 
bital elements that were assumed constant in the 
fit. On the other hand, projecting the orbital pa- 
rameters in Table 1 to those at a different epoch 
within the time span of the data set leads to much 
less variations in the orbital elements and appar- 
ent long term stability of the system (Marcy et al. 
2001). For some choices of epoch, the system is 
not in the 2:1 resonance, and it eventually becomes 
unstable (Laughlin & Chambers 2001). Clearly, a 
dynamical fit of system characteristics to the RV 
observations is necessary to constrain the param- 
eters that define an RV curve for GJ876. 

The publication of the data set by Marcy et 
al. (2001) allowed Laughlin & Chambers (2001) 
to perform such a fit. They assumed that the two 
planets are on coplanar orbits, and used two meth- 
ods to minimize xt as a function of the initial sys- 
tem parameters. Starting with the 2-Kepler fit, a 
Levenberg-Marquardt minimization scheme driv- 
ing an TV-body integrator was used to find a local 



minimum in xt- A second method uses a genetic 
algorithm combined with a simple model for the 
variations of the orbital elements due to the 2:1 or- 
bital resonance to search for the global minimum in 
xt- Both methods converged to similar solutions. 
The dynamical fits are sensitive to sini since the 
determined masses (and perturbations) grow as i 
is decreased. The best-fit solutions obtained using 
the Levenberg-Marquardt method are included in 
Table 1 for both the Keck+Lick data and Keck 
data alone. The Keck data is much more accu- 
rate than the Lick data and the solution for this 
data alone has a minimum xt — 1-59 with an rms 
scatter of 6.86 ms -1 for s'mi = 0.55. Although 
the Lick data is less accurate than the Keck data, 
the longer time span of observation that includes 
the Lick data may give the better solution. The 
solution for the Keck+Lick data has a minimum 
XI, = 1-46 with an rms scatter of 13.95 ms -1 for 
sin? = 0.78, but the xt minimum is broad, with 
several nearby solutions giving nearly as good a fit. 
The conclusion of Laughlin & Chambers is that 
the broad minimum in xt allows probable values 
of 0.5 < sin? < 0.8. (Rivera & Lissauer 2001 ex- 
amined mutually inclined orbits and got similarly 
small xt values for several different configurations, 
but they confirmed the Laughlin & Chambers re- 
sults for coplanar orbits.) 

Figure 1 shows the trajectories of the motions 
for 3100 days (the average periapse precession 
period for both planets) for the two Laughlin- 
Chambers best- fit solutions in plots of ej sin 8j ver- 
sus ej cos 9j , where ej is the eccentricity of the jth 
planet (with j = 1 and 2 for the inner and outer 
planets, respectively), 

01 = Ai — 2A2+TU1, , . 

2 =A 1 -2A 2 +tu 2 , 1 > 

are the two lowest order, eccentricity- type mean- 
motion resonance variables at the 2:1 commensu- 
rability, and Xj and Wj are the mean longitude 
and longitude of periapse of the jth planet. The 
two Laughlin-Chambers solutions place the sys- 
tem deep in both mean-motion resonances, with 
0i and 62 librating about 0° with remarkably small 
amplitudes. The simultaneous librations of B\ and 
62 about 0° mean that the secular resonance vari- 
able 

<9 3 = W! - W2 = 6>i - 2 (2) 
also librates about 0°. Although the parameters 
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Fig. 1. — Small amplitude librations of the two 2:1 
mean-motion resonance variables, 9\ = Ai — 2A 2 + 
w\ and 02 = Ai — 2A2+G72, about 0° for the GJ 876 
planets. Trajectories for 3100 days (the average 
periapse precession period for both planets) are 
shown in plots of ej sin 9j versus ej cos 0j for the 
Laughlin & Chambers (2001) best-fit solutions to 
the Keck data alone and the combined Keck and 
Lick data. 

will be more tightly defined as more data is ac- 
quired, the almost singular nature of the low am- 
plitude librations indicates that the real system 
is most likely indeed locked in multi-resonance li- 
brations at the 2:1 mean-motion commcnsurabil- 
ity. This means that the line of apsides of the 
two orbits are nearly aligned with conjunctions of 
the two planets always occurring very close to the 
periapse longitudes. The resonance configuration 
ensures that the planets can never make close ap- 
proaches in spite of their large masses, and barring 
external perturbations or an unreasonably high 
dissipation of tidal energy in the planets, the sys- 
tem should be stable for the lifetime of the star. In 
fact, we find that the librations are stable even for 
amplitudes of the inner planet resonance variable 
#i approaching 45° if we induce larger amplitudes 
of libration of the resonance variables by chang- 
ing the initial value of the mean anomaly of the 
inner planet. The existence of the mean- motion 
resonances means that the assumption that the 



Fig. 2. — Variations in the semi-major axes, a\ 
and d2, and eccentricities, e\ and e 2 , of the GJ 876 
planets for the Laughlin-Chambers Keck+Lick so- 
lution. The small amplitude librations of the res- 
onance variables ensure that the semi-major axes 
and eccentricities have little variation. 

orbits are nearly coplanar should be correct, and 
such coplanarity is consistent with the accretion of 
planetary bodies in a nebular disk surrounding the 
forming star. We shall assume that the orbits are 
coplanar hereinafter. Furthermore, there is not 
much to distinguish the two Laughlin-Chambers 
solutions in Figure 1, and we shall assume the pa- 
rameters appropriate to the solution based on both 
the Keck and Lick data. 

Figure 2 shows the variations in the semi-major 
axes and eccentricities of both planets for 10 4 days 
for the Laughlin-Chambers Keck+Lick solution. 
The small amplitude librations of the resonance 
variables about 0° ensure that both the eccentric- 
ities and semi-major axes have little variation in 
spite of the large mutual perturbations. This is 
why the 2-Kepler fit of Marcy et al. (2001) could 
produce a fit that was not too bad. 

The fact that both 2:1 mean-motion resonance 
variables librate about 0° in Figure 1 contrasts 
with the geometry of the Io-Europa 2:1 orbital 
resonances at Jupiter, where the resonance vari- 
able involving Io's longitude of periapse, 0\, li- 
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brates about 0° but the resonance variable in- 
volving Europa's longitude of periapse, 2 , librates 
about 180° (e.g., Peale 1999). This means that in 
the Io-Europa case the lines of apsides are anti- 
aligned, with the periapses 180° apart (i.e., 3 li- 
brates about 180°). Conjunctions occur when Io 
is near periapse and Europa is near apoapse. In §2 
we explain why the GJ 876 and Io-Europa systems 
have different resonance configurations by consid- 
ering a condition for stable simultaneous librations 
of the two mean-motion resonance variables. In 
particular, the longitudes of periapse should pre- 
cess in the same direction at the same average rate, 
so that the relative alignment of the lines of ap- 
sides is maintained. In the Io-Europa case, the 
orbital eccentricities are sufficiently small that the 
precession rate dzui/dt is dominated by a single 
term lowest order in eccentricities and containing 
cos 6^. The overall sign of the coefficient for this 
term is < for i = 1 (inner satellite) but > for 

1 = 2 (outer satellite), and retrograde precessions 
of both satellites require 0\ to librate about 0° and 

2 to librate about 180°. In the case of GJ876, 
the eccentricities are sufficiently large that there 
are large contributions to the precession rates from 
higher order terms whose cosine arguments are lin- 
ear combinations of the resonance variables. Sum- 
ming over all contributing terms, coincident ret- 
rograde precessions of the G J 876 planets require 
both 9i and 2 to librate about 0°. 

If the orbits of the two planets about GJ876 
were originally much further apart, with the ratio 
of their mean motions considerably greater than 
2:1, any process that drives the orbits toward each 
other could lead to capture into the 2:1 orbital 
resonances that we observe. In §3 we describe 
a particular migration process due to the gravi- 
tational interaction between the planets and the 
nebular disk surrounding the young GJ876. We 
note the simulations by Bryden et al. (2000) and 
Kley (2000), which show that two planets that are 
massive enough to open gaps individually in the 
gas disk can rather quickly clear out the disk ma- 
terial between them, if they are not separated too 
far. Disk material outside the outer planet ex- 
erts torques on the planet that are not opposed by 
disk material on the inside, and the outer planet 
migrates toward the star on the disk viscous time 
scale. Any disk material left on the inside of the in- 
ner planet exerts torques on the inner planet that 



push it away from the star. Thus the condition 
of approaching orbits necessary to form the res- 
onances is established. We discuss the effect of 
planet-nebula interaction on orbital eccentricities 
and note that the GJ876 system is in the inter- 
esting regime where it is uncertain whether eccen- 
tricity damping or growth is expected because of 
the rather large planet-star mass ratio of the outer 
planet. (Although the mass of the outer planet, 
M 2 , for the Laughlin-Chambers Keck+Lick solu- 
tion is only 2.4QM,,, M 2 /M = 7.17 x 10~ 3 because 
the stellar mass M = 0.32M Q .) 

In §4 we present the results of a series of nu- 
merical orbit integrations in which the orbits of 
the GJ 876 planets, initially far from the 2:1 mean- 
motion commcnsurability, are forced to approach 
each other. (The numerical methods are described 
in the Appendix.) We show that capture into 
the two 2:1 mean- motion resonances and the secu- 
lar resonance is certain if the orbital eccentricities 
start reasonably small and the rate of migration is 
not too fast. If there is no eccentricity damping, 
the eccentricities of both planets increase rapidly 
after resonance capture and exceed the observed 
values after a very short migration of the reso- 
nantly locked planets. Eccentricity damping of 
the form ej/e, oc dj/aj, where a dot over a sym- 
bol denotes d/dt and dj is the forced migration 
rate, leads to a termination in the eccentricity 
growth, and the eccentricities reach equilibrium 
values that remain constant for arbitrarily long 
migration in the resonances. We find that signif- 
icant eccentricity damping with |ei/e»| S> |dj/aj| 
is required to produce the observed eccentricities 
of the GJ876 system. In §5.1 we show that alter- 
native damping of the eccentricities by tidal dissi- 
pation within the star or planets is insignificant, 
and in §5.2 we discuss other related studies. Our 
conclusions are summarized in §6. 

2. COMPARISON WITH THE 
IO-EUROPA SYSTEM 

As shown in Figure 1, the GJ876 system has 
both lowest order, eccentricity-type mean-motion 
resonance variables at the 2:1 commcnsurability, 
0i = Ai — 2A2 + w\ and 2 = Ai — 2A2 + w 2 , 
and hence the secular resonance variable, #3 = 
zui — VJ2, librating about 0°. This resonance con- 
figuration is different from that of Io and Europa, 
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where 82 and 83 librate about 180°. It should 
be noted that the differences arc not due to the 
additional resonances involving Ganymede in the 
Io-Europa case. In the scenario where the reso- 
nances among the inner three Galilean satellites 
of Jupiter are assembled by differential tidal ex- 
pansion of the orbits (Yoder 1979; Yoder & Peale 
1981), Io is driven out most rapidly and the reso- 
nance variables involving Io and Europa only are 
captured into libration first. These resonance vari- 
ables have the same centers of libration before and 
after the 2:1 commensurability between Europa 
and Ganymede is encountered. The differences 
in the resonance configurations of GJ 876 and Io- 
Europa are instead due to the magnitudes of the 
eccentricities involved and can be understood from 
a condition for stable simultaneous librations of 
the two mean-motion resonance variables. 

This condition for stable simultaneous libra- 
tions of 61 and 82 is that the longitudes of pe- 
riapse, w\ and tzj 2 , on average precess at the same 
rate. For coplanar orbits, the equations for the 
variation of Wi, a,, and ej in Jacobi coordinates 
are (e.g., Peale 1986) 



dwi 
~~dT 

dai 
~~dt 

dfk 
dt 



dH 
fa~i dH 




(3) 
(4) 



dH 



Mleiy/iXidi dwi 



«r« (5) 



where M[ = M { YX=o M k/ Efe=o M ^ Mi = 

GM Mi/M-, M is the mass of the star (or 

Jupiter), Mi is the mass of the ith planet (or 
satellite), the Hamiltonian 



H = _j 2 GMoM i + ^ 

»=i 



2a; 



and the disturbing potential 

GM 1 M 2 „ 

$ = — - - GM M 2 

ri2 



1 

ro2 



1 

r-2 



(6) 
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If we neglect terms of order (Mi/M ) 2 and higher 
and assume that a x < a 2 , $ can be expanded to 



the form 



GM Y M 2 



a 2 
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(8) 



where 



»jkl 



= (J - fc)Ai - (J - £)X 2 + kw x - lw 2 , (9) 



and the summation is over the range < j < 00, 
—00 < k,£ < 00, and < m,n < 00. The coeffi- 
cient C 3 klmn {(5), where [3 — a\/a 2 , can be written 
in terms of the Laplace coefficient 6^ 2 (/3) and its 
derivatives (Brouwer & Clemence 1961). For a 
given cosine argument <pjke, the term lowest order 
in eccentricities is of order e[ fe 'e 2 £ '. 

Within the 2:1 mean- motion resonances and 
the secular resonance, the perturbations are dom- 
inated by the terms whose arguments are nearly 
fixed because of the resonances, i.e., those terms 
involving the resonance variables 61, 9 2 , #3, and 
their linear combinations. Since the terms in equa- 
tion (4) for dcii/dt and equation (5) for dei/dt 
are proportional to dH/dXi and dHjdwi and Aj 
and Wi appear in the cosine arguments only, the 
cosines are changed into sines and there is no sec- 
ular change in a, and the forced a if the resonance 
variables 9\, 9 2 , and 8 3 librate about either 0° or 
180°. (Note, however, that if the eccentricities are 
not small and dcii/dt and dei/dt are not dominated 
by a single term, there is also the possibility that 
the sums of all contributing terms are zero for res- 
onance variables having values other than 0° and 
180°.) 

We have calculated the contributions to the pre- 
cession rate dwi/dt from secular terms in equa- 
tion (8) with argument of the form kO^, (includ- 
ing the k = non-resonant secular term), up 
to fourth order in eccentricities, and from mean- 
motion resonance terms with argument of the form 
k9\ — £8 2 (k ^ £), up to third order in eccentric- 
ities, and they are listed in Table 2. We set the 
resonance variables to their average values, i.e., 
the libration center values. In addition, we ig- 
nore the small deviation of (3 from 2~ 2 / 3 (since 
the commensurability is not exact) and evaluate 
ClemnW at /? = 2- 2 /3. For Io and Europa, 
we adopt e\ — 0.0026 and e 2 = 0.0013, which 
are the equilibrium eccentricities before the res- 
onances with Ganymede are encountered in the 
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tidal scenario (Yoder & Peale 1981). The numeri- 
cal values in Table 2 are for the current semi-major 
axis of Io, but they can be scaled to other values 
of the semi-major axis since the precession rates 

3/2 

are proportional to a 1 ' ' . As we can see in Ta- 
ble 2, because the orbital eccentricities of Io and 
Europa are small, the precession rate dwi/dt is 
dominated by a single term lowest order in the 
eccentricities. This is the term with argument 9i 
and proportional to in the disturbing potential 
(eq. [8]), and whose contribution to dzui/dt is pro- 
portional to 1/a (eq. [3]). Since the coefficients of 
the ei cos 9\ and e 2 cos 9 2 terms are Ci 000 = —1.19 
and Cq_ 100 = +0.43, respectively, retrograde pre- 
cessions of both Io and Europa require 6\ = 0° and 
9 2 = 180°. Furthermore, the requirement that the 
regression rates are identical implies a simple re- 
lationship between the eccentricities: 



/Jni 



M 2 ClOOQ 

M ei 



. . Mi Co-ioo . • 

+ ^7 scc i — — n 2 — hra7 S cc,2, 

Mo e 2 



(10) 

where rij is the mean motion and the secu- 
lar motion vj scc .i includes contributions from 
secular terms in the disturbing potential and 
the additional secular motion induced by the 
oblatcness of Jupiter. As we can see in Ta- 
ble 2, the contributions from the secular terms 
in the disturbing potential and even the secu- 
lar motion induced by the oblateness of Jupiter 
(which is +0.12° day -1 for Io's orbit) are small 
compared with the <~ — 1.5°day _1 precession 
induced by the first order mean-motion reso- 
nance terms with arguments 9\ and 02- Thus 
ei/e 2 « -r 1/2 (C 1 2 000 /C 1 _ 100 )(M 2 /M 1 ) = 1.9, 
where the last equality is for the masses of Io and 
Europa. 

For GJ 876, we know from numerical orbit in- 
tegration of the Laughlin-Chambcrs Keck+Lick 
solution that on average dzui/dt = dw 2 /dt = 
—0.116° day -1 . The numerical values in Table 2 
are obtained using e\ — 0.255, e 2 = 0.035, and 
a\ = 0.130 AU, which are the average values for 
the Laughlin-Chambers Keck+Lick solution (see 
Fig. 2). As we would expect from the analysis 
above for the Io-Europa case, with 6\ = 9 2 =0°, 
the contributions to dwi/dt and dw 2 /dt from the 
d cos 9i and e 2 cos 9 2 terms, respectively, have op- 
posite signs. However, because the eccentricities 
are large, these terms no longer dominate the pre- 
cession rates. In particular, the largest contri- 



bution to dw 2 /dt comes from terms second or- 
der in eccentricities — the largest of which being 
the e\e 2 cos(#i + 9 2 ) term, which is of order ei(= 
0.255) lower than the e 2 cos 9 2 term, but whose 
coefficient is large and negative (C 3 _ 100 = —4.97). 
We can also see from Table 2 that both dw\/dt and 
dw 2 /dt converge only slowly with the order of the 
k9i — £9 2 terms. With the inclusion of terms up to 
e 3 , the precession rates have the correct, negative 
signs when 9\ = 9 2 = 9$ = 0°, but the magnitude 
of the precession rate for the inner (outer) planet 
is significantly larger (smaller) than that for both 
planets from numerical orbit integration. Based 
on the magnitudes and the alternating signs of the 
contributions from the lower order terms, the ex- 
pected contributions from the e 4 (and higher or- 
der) terms are consistent with their bringing the 
results into agreement with the actual precession 
rates. There is no simple analytic expression relat- 
ing ei to e 2 when the eccentricities are large, but 
as in the small eccentricity limit, the eccentricities 
are related by the requirement that the periapses 
precess at the same rate. 

Although the analysis in this section is able to 
explain why the GJ876 and Io-Europa systems 
have different resonance configurations, the slow 
convergence of the series for moderate to large ec- 
centricities means that it has limited usefulness 
if one is interested in the more general question 
of what stable configurations are possible for dif- 
ferent periapse precession rates and masses. A 
practical way to investigate this latter question 
is through numerical migration calculations like 
those in §4.1, which drive a system through a se- 
quence of configurations with different precession 
rates. As we shall see, for systems with masses 
like those in GJ876, there are stable configura- 
tions with #i, 9 2 , and 9$ librating about 0° for 
0-15 < ei < 0.86 (see Fig. 3). We have also 
found that configurations with resonance variables 
librating about angles other than 0° and 180° are 
possible when the masses are different from those 
in GJ876; these configurations will be discussed 
in a subsequent paper. 

3. MIGRATION SCENARIO FOR ORI- 
GIN OF RESONANCES 

We now turn to the question of the origin of the 
two 2:1 mean- motion resonances and the secular 
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resonance in the GJ 876 system. A scenario for 
the origin of this resonance configuration is that 
it was assembled by differential migration of or- 
bits that were initially much further apart. We 
shall see in §4 that libration of all three resonance 
variables is easily established by any process that 
drives the orbits toward each other. Although the 
results presented in §4 are quite general, they will 
be interpreted in the context of a particular dif- 
ferential migration process, namely via the grav- 
itational interaction between the planets and the 
nebula from which they form. 

ft is generally accepted that planets form in a 
disk of gas and dust that surrounds a young star. 
Since the two planets around GJ 876 are of order 
Jupiter mass and are presumably gas giants, they 
must have formed within the lifetime of the gas 
disk. If a planet is sufficiently massive, the torques 
exerted by the planet on the gas disk can open 
an annular gap in the disk about the planet's or- 
bit. The conditions for gap formation are that the 
Roche radius of the planet, tr = (M/3Af ) 1 ^ 3 a, 
exceeds the scale height H of the disk, or equiva- 
lently, 

M/M >3(ff/a) 3 = 3.75xl0- 4 ^) 3 , (11) 

and that the viscous condition 

M/M > AQv/{VLa 2 ) = A0a{H/af 

\A x 10- 3 / V0.05 J 

(12) 

is satisfied (e.g., Lin & Papaloizou 1993). In the 
above equations, M and M are the masses of the 
planet and the star, a is the semi-major axis of the 
planet's orbit, O is the angular Kepler speed at a, 
and the kinematic viscosity v is expressed using 
the Shakura-Sunyaev a prescription: v = aH 2 Q. 
Since the planet-star mass ratios of the planets 
around GJ 876 are M 1 /M = 2.28 x 1CT 3 and 
M 2 /M = 7.17 x 10 ~ 3 , these planets are expected 
to open gaps individually during their growth to 
their final masses if their orbits are sufficiently 
far apart, unless H/a > 0.09 or a 1 / 2 (H/a) > 
7.5 x 10~ 3 . 

The numerical values a = 4x 10~ 3 and H/a = 
0.05 are typical of models of protoplanetary disks 
(e.g., Bryden ct al. 2000; Kley 2000; Papaloizou 



et al. 2001). Hartmann et al. (1998) have also 
inferred from observed properties of T Tauri disks 
that a ~ 10~ 2 . The most promising source of an 
effective viscosity in accretion disks is magnetohy- 
drodynamic (MHD) turbulence initiated and sus- 
tained by the magnetorotational instability, which 
is capable of producing an effective a as large as 
0.1 (e.g., Stone et al. 2000). However, protoplan- 
etary disks are probably too weakly ionized for 
MHD turbulence to develop fully, except at small 
radii (< 0.1 AU) and possibly in a layer near the 
surface of the disk at larger radii (Gammie 1996). 
Another possible source of effective viscosity is 
damping of density waves excited by many small 
(of order Earth mass) planets, which is capable 
of producing an effective a < 10~ 3 (Goodman & 
Rafikov 2001). 

Bryden et al. (2000) and Kley (2000) have per- 
formed hydrodynamic simulations of a system con- 
sisting of a central star, a gas disk and two planets. 
The planets are massive enough (M/M ~ 10 -3 ) 
to open gaps individually, and the ratio of their ini- 
tial semi-major axes is a\/ai ~ 1/2. They found 
that the planets clear out nearly all of the nebular 
material between them in a few hundred orbital 
periods. Then the torques exerted by the nebu- 
lar material outside the outer planet drives that 
planet toward the star, whereas any nebular ma- 
terial left on the inside of the inner planet drives 
that planet away from the star. The depiction of 
the inner disk means that the inner planet may 
not move out very far, but the net effect is always 
to drive the orbits of two massive planets toward 
each other. The time scale on which the planets 
migrate is the disk viscous time scale, whose in- 
verse is (Ward 1997) 




(13) 

where a dot over a symbol denotes d/dt. 

In addition to its effect on the semi-major 
axis, planet-nebula interaction can also affect the 
orbital eccentricity of a planet (e.g., Goldreich 
& Tremaine 1980; Artymowicz 1992, 1993; Pa- 
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paloizou et al. 2001). A planet interacts with 
a disk in the vicinity of Lindblad and corotation 
resonances. The leading contribution to a is due 
to Lindblad resonances with the I = m Fourier 
components of the planet's perturbation potential. 
(The indices I and m for the Fourier series in time 
and azimuthal angle in this context should not be 
confused with I and m in eq. [8].) For a planet or- 
biting in a disk gap, the so-called co-orbital Lind- 
blad and corotation resonances are not important, 
and the leading contributions to e are damping 
due to I = m ± 1 corotation resonances and ex- 
citation due to I = m — 1 outer and I = m + 1 
inner Lindblad resonances. The net effect of the 
corotation resonance damping and the Lindblad 
resonance excitation depends on the distribution 
of the nebular material, which is itself determined 
by the interaction with the planet. If the gap is not 
too wide and many resonances of both types are 
present, previous calculations indicate that there 
is net eccentricity damping. For example, if we 
consider an outer disk of constant surface mass 
density £ that extends radially from a + A to oo, 
with A « a, then integration of equations (30) 
and (31) of Goldreich & Tremaine (1980) yields 




(Note that the above equations are identical to eqs. 
[109] and [110] of Goldreich & Tremaine where 
they set A = 2a/3ra max .) 

In the case of GJ876, since the planets clear 
out nearly all of the nebular material between 
them and the inner disk is likely to be depleted, 
the dominant resonant interactions are expected 
to be those between the outer disk and the outer 
planet. However, the estimate (14) for e is not 
adequate because, in addition to the fact that the 
mass distribution near the inner edge of the outer 
disk is not modeled properly, the outer planet is 
sufficiently massive that only a few low-m reso- 
nances are likely to be present in the disk. In 
fact, if the outer planet is able to open a gap out 
to the 2:1 commcnsurability, there would be no 
corotation resonances and only one Lindblad res- 
onance with I = m — 1 = 1 at the 3:1 commen- 
surability, and there would be eccentricity growth 



instead of damping (Artymowicz 1992). (In addi- 
tion, the disk can become eccentric and the growth 
of the planet's orbital eccentricity can be enhanced 
by the interaction with the eccentric disk, if the 
planetary mass is comparable to a characteristic 
mass of the disk; Papaloizou et al. 2001.) From 
a comparison of the resonant and viscous torques, 
an approximate condition for opening a gap out 
to the 2:1 commensurability is (Artymowicz 1992; 
Lin & Papaloizou 1993) 

M/Mq > 2.8a 1 / 2 {H/ a) 

„„ ln3 / a \ l/2 (H/a\ 
= 8.9 x 10" 3 - . 

(16) 

Since the planet-star mass ratio of the outer planet 
around GJ876 is M 2 /M = 7.17 x 10~ 3 and close 
to the critical value in equation (16), without 
knowing the exact values of the disk parameters, it 
is unclear whether one should expect eccentricity 
damping or growth. However, as we shall see in §4, 
significant eccentricity damping with |e/e| » |d/a| 
is required to produce the observed eccentricities 
of the GJ 876 system, unless the migration after 
resonance capture is severely limited. Therefore, 
at least the condition (16) must not be satisfied, 
implying that 

a 1/2 (H/a) > 0.36M 2 /Af 

, / M 2 /M \ 

= 2.6 x 10" 3 - — 

V7.17x 10- 3 J 

(17) 

for the outer disk of the young GJ 876 system. A 
more detailed analysis using hydrodynamic simu- 
lations is necessary to determine whether sufficient 
eccentricity damping can be produced by such an 
outer disk. 

4. NUMERICAL RESULTS FOR MI- 
GRATION SCENARIO 

In this section we present the results of a se- 
ries of numerical orbit integrations designed to de- 
termine the conditions under which the dynami- 
cal properties of the current GJ 876 system could 
be produced by any process (such as the planet- 
nebular interaction discussed in §3) that drives 
the orbits of the two planets toward each other. 
We consider a system consisting of a central star 
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and two planets, where the stellar and plane- 
tary masses are those for the Laughlin-Chambers 
Keck+Lick solution (Table 1). Unless stated oth- 
erwise, the planets are initially on coplanar, cir- 
cular orbits, with the mean longitudes differing 
by 180° and the ratio of the semi-major axes 
a\jai = 1/2, far from the 2:1 mean-motion com- 
mensurability. 

In addition to the mutual gravitational inter- 
actions of the star and the planets, we force the 
osculating semi-major axis of planet i to migrate 
at a rate dj. In most cases, we assume that only 
the outer planet is forced to migrate inward, with 
a migration rate of the form d 2 /a 2 = constant. 
The effects of adopting a more general form of 
the migration rate (e.g., 0,2/0,2 being a function 
of a 2 ) or forcing the inner planet to migrate out- 
ward are discussed in §4.2. For the calculations in 
§4.2, we also damp the osculating eccentricity at 
a rate ej. We adopt a damping rate of the form 
= —K\ai/ai\, where if is a positive constant. 
As we shall see, this form of eccentricity damping 
has the convenient property that the eccentricities 
reach equilibrium values after capture into the 2:1 
resonances. We note that the relation e/eocd/ais 
valid for the simple estimates (14) and (15) for the 
interaction of a massive planet with an outer disk 
if A/a w constant and that the same relation be- 
tween the damping and migration rates also holds 
for a planet that is too small to open a gap in a 
gas disk with constant H/ a and undergoes the so- 
called type I migration (Artymowicz 1993; Ward 
1997). ' 

The numerical orbit integrations were per- 
formed using the symplectic integrator SyMBA 
(Duncan, Levison, & Lee 1998) modified to in- 
clude the orbital migration and eccentricity damp- 
ing terms. SyMBA is based on a variant of the 
Wisdom-Holman (1991) method and employs a 
multiple time step technique to handle close en- 
counters. (The latter feature is not essential for 
the integrations presented here.) We confirmed 
the results by repeating some of the integrations 
using a Bulirsch-Stoer integrator that has also 
been modified to include migration and damping. 
We also checked that both modified integrators 
give the correct exponential decays in a and e 
when they are used to integrate the orbit of a 
single planet with a/ a = constant and e/e = 
constant. We describe how the integrators were 



modified in the Appendix. 

4.1. Migration without Eccentricity Damp- 
ing 

We consider first inward migration of the outer 
planet without eccentricity damping. In Figure 3 
we show the results of a calculation with a\ = 
0.5 AU and 02 = 1.0 AU initially and 0,2/02 = 
—5 x 10 _5 yr _1 . The migration rate is consis- 
tent with that due to planet-nebular interaction 
at ~ 1AU (cq. [13]). Figure 3a shows the time 
evolution of the semi-major axes and eccentrici- 
ties, and Figure 36 shows the evolution of the two 
2:1 mean-motion resonance variables, 9\ and 02, 
in plots of eisin^ versus ejcosflj. Initially, only 
the outer planet migrates inward at the prescribed 
rate. When the 2:1 mean-motion commensurabil- 
ity is encountered, both 6\ and 62 (and hence the 
secular resonance variable #3) are captured into li- 
bration about 0°. We do not see libration of 82 
and 63 about 180°, which is expected for small ec- 
centricities (see §2), because the planetary masses 
are so large that fairly large eccentricities (with 
d reaching about 0.15) are generated before the 
2:1 commensurability is encountered. Because of 
the forced migration, the centers of libration are 
actually slightly offset from 0°. The resulting res- 
onant interaction slows down the migration of the 
outer planet and forces the inner planet to mi- 
grate inward while keeping the ratio of the semi- 
major axes nearly constant (01/02 ~ 2~ 2 / 3 ); it 
also causes the eccentricities to increase rapidly. 
The centers of libration remain near 0° and the 
amplitudes of libration remain small as the eccen- 
tricities increase. 

The sequence of configurations with increasing 
eccentricities that the system is driven through af- 
ter resonance capture is in fact a sequence with in- 
creasingly less negative periapse precession rates. 
As we discussed in §2, the forced eccentricities 
(and the libration centers) for a system with sta- 
ble simultaneous librations of 0\ and 9 2 are de- 
termined by the requirement that the longitudes 
of periapse on average precess at the same rate 
(which in turn is determined by Ai — 2A 2 + Wi — 
0). Although a longer integration indicates that 
the system does eventually become unstable when 
ei w 0.86, the existence of stably librating con- 
figurations with ei up to 0.86 is remarkable. In 
particular, the configurations with e\ > 0.71 have 
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a/a, = a/a 2 = -5xl(T 5 yr 1 

e/e, = e 2 /e 2 = 




2xl0 4 4xl0 4 



t (yr) 



prograde periapse precessions. We have integrated 
the configurations at t = 2 x 10 4 yr (with e\ w 0.59 
and retrograde periapse precessions) and 5 x 10 yr 
(with ei w 0.81 and prograde periapse preces- 
sions) forward with migration turned off and con- 
firmed that the system remains stable as seemingly 
indefinitely repeating configurations. 

As we can see in Figure 3 a, without eccentric- 
ity damping, the eccentricities exceed the observed 
values for the GJ 876 planets (dashed lines) when 
the semi-major axes of the resonantly locked plan- 
ets have decreased by only w 7% after capture into 
the resonances. This result is insensitive to the 
adopted parameters. If the migration rate is not 
too fast (see below) , the evolution of a system with 
different initial a 2 (but for convenience, the same 
0,1/0,2) or 0,2/02 is essentially identical to that 
shown in Figure 3 if we plot the semi-major axes in 
units of initial 02 and time in units of (0,2/0,2) _1 - 
Therefore, unless by coincidence resonance cap- 
ture occurs just before migration stops because of, 
e.g., nebula dispersal, eccentricity damping is nec- 
essary to produce the observed eccentricities of the 
G J 876 system. 

From a set of calculations similar to that shown 
in Figure 3, with initial ai/02 = 1/2, 02 = 1, 2, 
and 4AU, and different a 2 /a 2 , we find that cer- 
tain capture of both 2:1 mean-motion resonance 
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Fig. 3. — (a) Time evolution of the semi-major 
axes and eccentricities and (b) evolution of the 
mean-motion resonance variables 6\ = Ai — 2A 2 + 
mi and 62 = Ai — 2A 2 + W2 in plots of e,sin#j 
versus e» cos 9i for a calculation where the outer 
planet is forced to migrate inward with a 2 /a 2 = 
—5 x 10 -5 yr -1 and there is no eccentricity damp- 
ing. Both 9i and 62 are captured into small am- 
plitude libration about 0°, but the eccentricities 
exceed the observed values for the GJ876 plan- 
ets (dashed lines) shortly after resonance capture, 
when the semi-major axes of the resonantly locked 
planets decrease by only 7%. 

variables requires 

where a 2 is the semi-major axis of the outer planet 
when the 2:1 commensurability is encountered. 
For migration rate within a factor of a few of, 
but below, the above limit, both resonance vari- 
ables are captured into libration, but the centers 
of libration can be significantly different from 0° 
and the amplitudes of libration can be large. The 
condition (18) is satisfied by a migration rate due 
to planet-nebular interaction with the nominal pa- 
rameter values in equation (13) by almost 3 orders 
of magnitude. 

To examine the effects of orbital eccentricities 
on capture into the 2:1 resonances, we performed 
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a series of calculations with non-zero initial ec- 
centricities. Four types of initial conditions were 
considered: the outer planet was initially at cither 
periapse or apoapse, and either the initial e\ or e 2 
was fixed at 0.01 and the other initial eccentricity 
was varied. The inner planet was always started 
at periapse, with its mean longitude differing from 
that of the outer planet by 180°, and the remain- 
ing parameters were identical to those used in the 
calculation shown in Figure 3. Since gravitational 
interaction between the planets causes the eccen- 
tricities to fluctuate even when a\ and a 2 are close 
to their initial values, the eccentricities quoted be- 
low are the maximum eccentricities when a\ and 
a 2 are close to their initial values and not the ini- 
tial eccentricities. We find that certain capture of 
both 2:1 mean-motion resonance variables requires 

d < 0.06 and e 2 < 0.03. (19) 

For eccentricities above these limits, there is non- 
zero probability for the planets to be captured into 
higher order resonances (e.g., 5:2) encountered be- 
fore the 2:1 commensurability. 

4.2. Migration with Eccentricity Damping 

Figure 4 shows the results of a calculation sim- 
ilar to that shown in Figure 3, but with eccen- 
tricity damping of the form e 2 /e 2 = — K\a 2 /a 2 \, 
where K = 100. The capture of the resonance 
variables Q\ and 9 2 into libration about 0° and the 



Fig. 4. — Same as Fig. 3, but for a calculation 
with eccentricity damping of the form e 2 /e 2 = 
— K\a 2 /a 2 \, where K = 100. After resonance 
capture, in addition to librations of Q\ and 9 2 
about 0°, the eccentricities reach equilibrium val- 
ues close to the observed values for the G J 876 
planets {dashed lines) and remain constant for ar- 
bitrarily long migration in the resonances. (The 
jagged nature of the plot (b) prior to equilibrium 
is due to sparse sampling.) 

initial evolution after resonance capture are simi- 
lar to the case without damping. However, the ec- 
centricity growth eventually terminates when the 
damping balances the excitation due to resonant 
interaction between the planets. The eccentricities 
reach equilibrium values that remain constant for 
arbitrarily long migration in the resonances. With 
K = 100, the equilibrium eccentricities are close to 
the observed eccentricities of the GJ 876 system. 
At the end of the calculation shown in Figure 4, 
when t — 4.6 x 10 4 yr, the semi-major axes are also 
similar to those of the GJ 876 planets. 

We have integrated the configuration at the end 
of the calculation shown in Figure 4 forward with 
migration and damping turned off. The system re- 
mains stable, with almost no change in the ampli- 
tudes of libration (compare Figs. 46 and 5), which 
are somewhat smaller than those of the G J 876 
system (compare Figs. 1 and 5). Although the 
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Fig. 5. — Continued small amplitude librations of 
the mean-motion resonance variables 9\ and #2 af- 
ter termination of planet migration and eccentric- 
ity damping. The configuration at the end of the 
calculation shown in Fig. 4 is integrated forward 
for 500 yr with migration and damping turned off. 



libration amplitudes are fairly similar in the two 
Laughlin-Chambers best-fit solutions, we may find 
in the future, as more data is acquired, that im- 
proved best-fit solutions will have smaller libra- 
tion amplitudes. Alternatively, larger amplitudes 
could be generated after termination of migration 
and damping by encounters with remaining plan- 
ctcsimals. 

In Figure 6 we show the time evolutions of the 
eccentricities for a set of calculations with 02/02 = 
—5 x 10~ 5 yr -1 , e 2 /e 2 = — K\ci2/a2\, and different 
K. The equilibrium eccentricities decrease with 
increasing K and are significantly different from 
the observed eccentricities of the G J 876 system if 
K is more than a factor 2-3 larger or smaller than 
K = 100. 

As long as e2/e 2 = — K\ib2/a2\, where K is 
a positive constant, the equilibrium eccentricities 
are determined by the value of K only and are 
insensitive to either the magnitude or functional 
form of 0,2/0,2- This is demonstrated in Figure 7, 
where we plot the results of a calculation with the 



Fig. 6. — Decrease of the equilibrium eccentricities 
with increasing eccentricity damping rate. Time 
evolutions of the eccentricities for a set of calcu- 
lations with 02/02 — —5 x 10~ 5 yr -1 , e2/e2 = 
— K\a2/o2\, and different K are shown. 

same K{= 100) as that shown in Figure 4 but 
with a different form of the migration rate: a 2 = 
constant or a 2 /a 2 oc 1 /a 2 . The equilibrium eccen- 
tricities are identical in Figures 4 and 7. If &2 / &2 is 
not exactly proportional to 1 02/02!, the eccentrici- 
ties would decrease slowly (after reaching maxima) 
or increase slowly as the resonantly locked planets 
migrate, but the ratio of damping to migration 
rates, |(e 2 /e2)/(a2/a2)|, just before migration and 
damping stop should be close to 100 for the final 
eccentricities to be close to the observed values of 
the GJ 876 system. 

Thus far we have considered forced inward mi- 
gration of the outer planet only. This is the most 
likely situation in the planet-nebular interaction 
scenario discussed in §3, since the inner disk is 
likely to be depleted and the dominant interac- 
tions are expected to be those between the outer 
disk and the outer planet. However, if the inner 
disk is not depleted, at least initially, or migration 
and damping are due to another process, the in- 
ner planet may also be forced to migrate outward. 
In Figure 8 we show the results of a calculation 
where, for simplicity, the inner planet is forced 
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e 1 /e 1 = e 2 /e 2 = -100|a 2 /a 2 | 
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Fig. 7. — An example demonstrating that the 
equilibrium eccentricities are determined by the 
ratio of damping to migration rates, K, and do 
not depend on the functional form of the migra- 
tion rate. The results of a calculation with the 
same ratio (K = 100) as that shown in Fig. 4 but 
with a different form of the migration rate (a,2 = 
constant) are shown. 

to migrate outward at the same rate that the 
outer planet is forced to migrate inward (di ja\ = 
— d 2 /a 2 = 5x 10~ 5 yr _1 ) and ej/ej = —K\di/ai\. 
Initially, the inner (outer) planet migrates outward 
(inward) at the prescribed rate. After resonance 
capture, resonant interaction overcomes the forced 
outward migration of the inner planet, and both 
planets migrate inward slowly. The equilibrium 
eccentricities are close to the observed eccentrici- 
ties of the GJ 876 system when K = 10. Therefore, 
even with eccentricity damping of both planets, 
significant eccentricity damping with | ej/ej | 3> 
|dj/oi| is required to produce the observed orbital 
eccentricities of the GJ 876 planets. 

5. DISCUSSION 

5.1. Eccentricity Damping by Tidal Dissi- 
pation in the Star and Planets 

We have found in §4 that significant eccentricity 
damping with | e» /e» | 3> |dj/eij| is required to pro- 
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Fig. 8. — Time evolution of the semi-major axes 
and eccentricities for a calculation where the in- 
ner planet is also forced to migrate outward, 
with di/ai = —02/02 = 5 x 10 _5 yr _1 , ej/ej = 
— K\di/ai\, and K = 10. Even with eccentricity 
damping of both planets, K 10 is required to 
produce the observed eccentricities of the GJ 876 
planets {dashed lines). 

duce the observed eccentricities of the G J 876 sys- 
tem if the migration has been at all extensive after 
resonance capture. As we discussed in §3, it is as 
yet unclear whether sufficient damping could be 
produced by planet-nebula interaction, even if the 
condition (17) is satisfied and eccentricity damp- 
ing is expected (see also §5.2). In this subsection 
we show that alternative eccentricity damping by 
tidal dissipation within the star and planets dur- 
ing planet migration is completely negligible. 

The star GJ 876 would most likely still be in 
its pre-main-sequence contracting phase during 
disk evolution and planet migration. To esti- 
mate the rate of circularization of an orbit due 
to tidal dissipation in the star, we use the stel- 
lar radius when the time t since initial contrac- 
tion is approximately the migration time scale of 
\a/a\ — 2x 10 4 yr (eq. [13]) and keep the planets at 
their current distances from the star. The planets 
would have migrated for longer than the migration 
time scale (see, e.g., Fig. 4) and it would take time 
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for the planets to form. But by adopting this min- 
imum time to yield the maximum probable stellar 
radius and by keeping the planets at their current 
distances, we maximize our estimate of the tidal 
rate of circularization of an orbit. 

Since the mass of GJ876 is only 0.32M Q , we 
assume that GJ876 follows the nearly vertical 
Hayashi track in the HR diagram and remains fully 
convective during its contracting phase. Hence the 
star remains a polytrope of index n = 1.5, with ef- 
fective temperature T e w 3500 K (e.g., D'Antona 
& Mazzitelli 1994). From the virial theorem, the 
stellar luminosity 
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where a is the Stefan-Boltzmann constant, G is 
the gravitational constant, Mo and Ro are the stel- 
lar mass and radius, and the form of the gravita- 
tional energy for a polytrope of index n is given 
by Chandrasekhar (1939). Integration of equation 
(20), with 3/2(5 - n) = 3/7 for a polytrope of 
index n = 1.5, gives Ro s=s 3.85 x 10 11 cm after 
t = 2 x 10 4 yr, leading to a ratio of the stellar ra- 
dius to the orbital semi-major axis Ro/a w 0.197 
for the inner planet and Ro/a w 0.125 for the outer 
planet. 

Since the present orbital period of the inner 
planet of the GJ876 system is about 30 days, 
which is long compared to the periods of free oscil- 
lation of the star, we can use the rate for circular- 
izing an orbit corresponding to dissipation domi- 
nated by the equilibrium tide for the star, which 
is given by (Zahn 1989) 
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where q = M/M is the planet-star mass ratio 
and t f = (MoRq/L) 1 / 3 . We neglect the reduc- 
tion of the turbulent viscosity when the tidal pe- 
riod is short compared to the turnover time of the 
largest convective eddies and adopt the maximum 
A c ; rc 0.048 for a fully convective star. Substitu- 
tion of Ro/a from the previous paragraph, along 
with q for the Laughlin-Chambers Keck+Lick so- 
lution (Table 1), t f = (Af /47rCTT ( 4 ) 1 /3 ~ 0.574yr, 



and A c i rc w 0.048, into equation (21) gives 
w 9.1 x 10 _9 yr _1 for inner planet, 

trirc 



w 7.6 x 10 10 yr 1 for outer planet. 

(22) 

Since we have assumed that the planets are at 
their closest proximity to the star for the whole 
time the tidal effects are damping the eccentric- 
ities and that the star is inflated to a maximum 
size, these values are extreme upper bounds on 
the rate of eccentricity damping by tides raised on 
the star. Even in its inflated state, the dissipation 
in the star can have essentially no effect on the 
orbital eccentricities. 

The theoretical circularization rate in equation 

(21) has been tested by comparing with the ob- 
served circularization rates of binaries. It is in 
good agreement with the observed rates for bi- 
naries containing giant stars (Verbunt & Phinney 
1995), but even with the maximum A c i rc , it is 
about 50-100 times slower than the observed rates 
for binaries containing main-sequence solar-type 
stars, which have radiative cores (Claret & Cunha 
1997; Goodman & Oh 1997). It is not yet clear 
why there is a discrepancy in the latter case or 
that this discrepancy is relevant for a fully convec- 
tive pre-main-sequence star. Nevertheless, even if 
we increase the circularization rates in equation 

(22) by a factor 100, they are still much smaller 
than the migration rate due to planet-nebula in- 
teraction. 

Tidal dissipation within a gaseous planet damps 
its orbital eccentricity at a rate given by (e.g., 
Peale, Cassen, & Reynolds 1980) 
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where R, ki, and Q are the radius, the poten- 
tial Love number, and the dissipation function of 
the planet. If we adopt values for R, fc 2 , and Q 
similar to those of Jupiter, with R — 7 x 10 4 km, 
k 2 = 0.38 (Gavrilov & Zharkov 1977), and Q = 
5 x 10 4 (which is approximately the lower bound 
on Q for Jupiter; Yoder & Peale 1981), |e/e| w 
1.5 x 10" 12 y^ 1 for the inner planet of the GJ 876 
system. The rate for damping the outer planet's 
orbital eccentricity from dissipation within itself 
is of course even smaller. The planets would most 
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likely be larger and contracting during their mi- 
gration. However, unless the planetary radii are 
unrealistically large and comparable to the Roche 
radii, the eccentricity damping rate due to tidal 
dissipation within the planets is smaller than the 
migration rate due to planet-nebula interaction by 
several orders of magnitude. 

5.2. Other Studies 

After completing our calculations (Lee & Peale 
2001), two papers with complementary calcu- 
lations came to our attention (Sncllgrove, Pa- 
paloizou, & Nelson 2001; Murray, Paskowitz, & 
Holman 2001). 

Snellgrove et al. (2001) also find from numer- 
ical orbit integrations with forced migration and 
eccentricity damping of the outer planet that the 
orbital eccentricities of the GJ 876 planets require 
a short time scale for eccentricity damping com- 
pared to the migration time scale. Note, however, 
that they adopt the minimum planetary masses 
from the two-Kepler fit by Marcy et al. (2001) 
and, not surprisingly, have difficulties matching 
both ei and e 2 from their calculations to those 
from the same fit, since the eccentricities from 
the two-Kepler fit are not consistent with small- 
amplitude simultaneous librations of 8\ and 82- 
Sncllgrove et al. also develop an analytic reso- 
nance theory that is first order in the eccentrici- 
ties. As we showed in §2, a first order theory is 
inadequate for understanding the current GJ876 
resonance configuration, which has 0±, 82, and #3 
all librating about 0°. Furthermore, as we found in 
§4 (see, e.g., Fig. 3), because the eccentricities are 
excited to large values in the GJ 876 evolution be- 
fore resonance capture, there is no time during the 
evolution when a first order theory can be useful 
for this system. Thus all predictions coming from 
this theory should be viewed with caution. Snell- 
grove et al. also present a hydrodynamic simula- 
tion of the planet-nebula interaction, where both 
planets are inside a cavity with almost no disk ma- 
terial. At the end of this simulation, e\ w 0.34 and 
the semi-major axes of the planets have decreased 
by about 13% since resonance capture. We have 
performed a numerical orbit integration without 
eccentricity damping similar to that shown in Fig- 
ure 3 but with the planetary masses adopted by 
Sncllgrove et al. and find that e\ rts 0.38 for the 
same reduction in semi-major axes. Thus it is not 



clear that the eccentricities have reached equilib- 
rium values at the end of this simulation; even if 
they have, these equilibrium values are too large 
for the GJ 876 system. It appears that the disk 
model used in this hydrodynamic simulation is not 
very effective in damping the eccentricities, which 
emphasizes the uncertainty in such damping. 

Murray et al. (2001) consider mainly an al- 
ternative scenario in which the migration and ec- 
centricity damping of the outer planet are due to 
scattering of planetesimals in the disk population. 
Their numerical simulations confirm our results 
(and those of Sncllgrove et al. 2001) of easy cap- 
ture of an inner planet into resonance, the dual mi- 
gration of both planets in the resonance, and the 
growth of the eccentricities, although they limit 
their numerical and analytical studies to an outer 
planet of Jupiter mass and inner planets of several 
Earth masses and devote much of the discussion 
to resonances of higher order than those at the 
2:1 mean- motion commcnsurability. Thus their 
results are not directly applicable to the GJ876 
system. 

6. CONCLUSIONS 

Radial velocity measurements by Marcy et al. 
(2001) have revealed two planets in resonant or- 
bits about the star GJ876. The remarkable or- 
bital fit obtained by Laughlin & Chambers (2001), 
which finds both lowest order, eccentricity-type 
mean- motion resonance variables at the 2:1 com- 
mensurability librating with small amplitudes, 
means that the resonances are almost certainly 
real and indefinitely stable. The existence of the 
eccentricity-type resonances implies that the as- 
sumed coplanarity of the orbits is probably close 
to reality. 

The GJ876 planetary system has revealed 
properties of the 2:1 orbital resonances that have 
not been observed nor analyzed before. The li- 
bration of both lowest order mean-motion res- 
onance variables, 6\ = Xi — 2X2 + wi and 
82 = Ai — 2X2 + ZU2, and the secular resonance 
variable, 83 = w\ — W2, about 0° was not antici- 
pated, as the familiar Io-Europa 2:1 resonance has 
8\ librating about 0°, but 82 and 83 librating about 
180° — a configuration that would persist in the 
absence of Ganymede. Thus conjunctions for the 
Jovian satellites occur when Io is near periapsc 
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and Europa is near apoapse, whereas conjunc- 
tions of the two planets about G J 876 occur when 
both planets are near periapse. We understood 
this to be mainly a function of the eccentricities of 
the orbits, where the resonance configuration with 
6»i w 0° and 6 2 ~ 3 ~ 180° must obtain when the 
eccentricities are small, but the resonance configu- 
ration with 9\ w 9 2 ~ #3 ~ 0° prevails for a system 
with masses like those in GJ 876 when the eccen- 
tricities are large. A necessary condition for stable 
simultaneous librations of both mean-motion res- 
onance variables is that w\ — w 2 on average, so 
that the relative alignment of the lines of apsides 
of the two orbits is maintained. The periapse pre- 
cessions are dominated by resonant terms in the 
disturbing potential whose arguments arc 9\, 9 2l 
#3, and their linear combinations. The dominance 
of the lowest order terms when the eccentrici- 
ties are small allows equal precession rates only 
if the lines of apsides are anti-aligned (61 w 0°, 
9 2 rts #3 s=a 180°), whereas dominance of higher or- 
der terms when the eccentricities are large results 
in equal precession rates for a system with masses 
like those in GJ 876 when 9 1 w 9 2 w 63 « 0°. The 
equality of the precession rates also determines a 
relationship between the eccentricities of the two 
orbits, although there is no simple analytic expres- 
sion for this relationship when the eccentricities 
are large, since the first order theory is not a good 
representation. The stability of the GJ876 reso- 
nance configuration for values of e\ up to 0.86 was 
also a surprise. 

Any process that drives the two originally 
widely separated orbits toward each other can 
result in capture of the planets into orbital reso- 
nances at the 2:1 commensurability. The naturally 
occurring situation where nebular disk material is 
cleared between two planets sufficiently massive 
to individually open gaps in the disk (Bryden 
ct al. 2000; Kley 2000) leads to inward migra- 
tion of the outer planet and possibly outward 
migration of the inner planet. Thus the likely 
origin of the resonances in the GJ 876 system is 
this differential planet migration due to torques 
induced by the planet-nebula interaction. We 
have shown that forced inward migration of the 
outer planet of the GJ 876 system results in cer- 
tain capture of 61, 9 2 and hence 9 3 into libra- 
tion if initially e± < 0.06 and e 2 < 0.03 and 
1 62/02 1 < 3 x 10- 2 (a 2 /AU)- 3 / 2 yr- 1 . The lat- 



ter rate is three orders of magnitude higher than 
the likely rate of - 5 x 10- 5 (a 2 / AU)" 3 / 2 yr" 1 
due to planet-nebular interaction. The bounds 
on the eccentricities result not so much from the 
transition from certain to probabilistic capture at 
the 2:1 resonances but from likely capture into 
higher order resonances such as 5:2 before the 2:1 
commensurability is encountered. 

Continued migration of the planets while locked 
in the 2:1 resonances leads to rapid growth in the 
orbital eccentricities that exceed the observed ec- 
centricities of the G J 876 system after only a fur- 
ther decrease in the semi-major axes of about 7% 
if there is no eccentricity damping. So unless 
resonance capture occurred near the end of mi- 
gration, the observed values of the eccentricities 
require eccentricity damping. With damping of 
the form ej/e, = —K\a,i/a,i\, where if is a posi- 
tive constant, eccentricity growth is terminated at 
values of the eccentricities that increase with de- 
creasing K, and the eccentricities remain constant 
for indefinite duration of the migration. The ob- 
served eccentricities result for K w 100 if there is 
forced migration and eccentricity damping of the 
outer planet only, but for K « 10 if there is also 
forced migration and eccentricity damping of the 
inner planet. This result is independent of the 
magnitude or functional form of hi/cii as long as 
ej/e, = —K\hi/ai\ is preserved. Relaxing the last 
condition leads to a slow drift in the eccentricities 
during migration and would require the migration 
to terminate as the eccentricities pass through the 
observed values. 

Existing analytic estimates of the effects of 
planet-nebular interaction are consistent with ec- 
centricity damping of the form e^/e^ = —K\hi/ai\, 
if the planet-star mass ratio is not too large (e.g., 
Goldreich & Tremaine 1980; Artymowicz 1992, 
1993; Lin & Papaloizou 1993; Ward 1997). How- 
ever, the planet-star mass ratio of the outer planet 
of the GJ 876 system is sufficiently close to the 
critical value separating eccentricity growth from 
damping for nominal values of the disk param- 
eters that it is uncertain whether such damping 
would occur. Even if the disk parameters are 
such that eccentricity damping would occur, it 
is not clear that the magnitude of K would be 
sufficiently large to constrain the eccentricities in 
the GJ876 system to the observed values. We 
have shown that the alternative eccentricity damp- 
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ing mechanism involving the dissipation of tidal 
energy within the star and the planets is com- 
pletely negligible. Further long-term hydrody- 
namic simulations with different physical assump- 
tions and parameters are required to determine 
whether planet-nebular interaction could produce 
sufficient eccentricity damping to allow arbitrary 
migration of the planets within the resonances in 
the young GJ 876 system while preserving eccen- 
tricities comparable to those observed. If not, the 
migration must have been finely tuned to stop 
when the system had progressed to its observed 
state, although this latter constraint is too ad hoc 
to be believable. 

It is a pleasure to thank Greg Laughlin for fur- 
nishing the details of the dynamical fit by Laughlin 
& Chambers before publication and Lars Bildsten 
for pointing out the enhancement of tidal dissipa- 
tion in a pre-main-sequence star. We also thank 
D.N.C. Lin, J.J. Lissaucr, and W.R. Ward for 
informative discussions. This research was sup- 
ported in part by NASA grants PGG NAG5-3646 
and OSS NAG5-7177. 
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A. NUMERICAL METHODS 



In this appendix we describe how the symplectic integrator SyMBA and the Bulirsch-Stoer integrator used 
for the numerical orbit integrations presented in §4 were modified to include the forced orbital migration 
and eccentricity damping terms. 

A second-order symplectic integrator for a Hamiltonian of the form H = H + Hi, where H and Hi are 
separately integrable, can be represented as 

E QEi(t)E Q. (Al) 

The three operators in equation (Al) represent a single step of an algorithm that starts with evolving the 
system under the influence of H only for half a time step r/2, then evolving it for a full time step r 
under the influence of Hi, and then evolving it for another half a time step r/2 under Hq. For example, 
in the Wisdom-Holman (1991) method for the gravitational A-body problem, the Hamiltonian in Jacobi 
coordinates is divided into H that describes the Keplerian motion of the planets around a central star 
and Hi that describes the perturbation of the planets on one another and is a function of the positions 
q i only. A step of the Wisdom-Holman method is thus: (1) Each planet evolves along a Kepler orbit for 
time r/2; (2) each planet receives a kick to its momentum of the amount —TdH\jdq i while q i is unchanged 
(since Hi does not involve the canonical momenta); (3) each planet evolves along a Kepler orbit for time 
r/2, starting with the new momentum after the kick. Recursive application of the basic algorithm (Al) 
allows one to construct symplectic integrators for Hamiltonians that consist of more than two integrable 
parts. For example, a single step of an algorithm for a Hamiltonian of the form H = Ho + Hi + H2 is 
E q (t/2)Ei(t/2)E 2 (t)Ei(t/2)E q (t/2). 

The symplectic integrator SyMBA (Duncan ct al. 1998) is based on a variant of the Wisdom-Holman 
method, with the gravitational A-body Hamiltonian written in terms of positions relative to a central star 
and barycentric momenta, and employs a multiple time step technique to handle close encounters. In the 
SyMBA algorithm, the Hamiltonian is divided into more than two parts and the recursive application of 
the algorithm (Al) discussed in the previous paragraph is utilized. Although the additional forced orbital 
migration and eccentricity damping terms are not Hamiltonian, they can be incorporated in an analogous 
manner. Our modified algorithm is 

Ea (£) E e (I) £ grav (r) E e (I) E a (I) , (A2) 

where E glav (T) denotes a complete time step for the conservative gravitational TV-body problem using the 
SyMBA algorithm, and E a {r/2) and E e (r/2) denote changing the canonical variables according to the 
imposed d, and terms, respectively, for time r/2. During the application of the d, term, all of the other 
orbital elements are constant, and a^i = a^o exp(rdj/2aj), where a^o and a^i are <Xj at the beginning and 
end of the step, if di/a, = constant (this can be easily generalized for, e.g., dj/aj oc a]). Note that we do 
not use truncated approximation such as a^i = a^ (l +rd i /2a i ). Similarly, the E e (r/2) step changes the 
eccentricities according to e^i = ej ; o exp(— rK\di/ai\/2) if e^/ej = — K\di/ai\. By modifying the algorithm 
in a symmetric manner and using exact solutions in the E a (r/2) and E e (r/2) parts, there should be little 
(if any) secular growth in the energy error (e.g., Mikkola 1998). 

For the Bulirsch-Stoer integrator, additional terms must be included in the equations of motion in Carte- 
sian coordinates to account for the forced orbital migration and eccentricity damping. In the following, 
we simplify the notation by considering a specific planet and dropping the subscript. Let (x,y, z) be the 
Cartesian coordinates of the planet with respect to the star and r be the distance of the planet from the 
star. The osculating orbital elements a, e, i, f, ui, and f2 are the semi-major axis, eccentricity, inclination, 
true anomaly, argument of periapse, and longitude of the ascending node on the xy plane, respectively. The 
additional terms in the equations of motion due to the forced migration d and eccentricity damping e terms 
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arc: 



dx 

~dl 

dx 
~dt 



dx 
~dt 

dx 
~dt 



dx . dx 



da 
dx 



de 



dx 



da U de 6 ' 



(A3) 



(A4) 



with similar expressions for the other coordinates. To evaluate the partial derivatives in equations (A3) and 
(A4) and similar expressions for the other coordinates, we need to express the position and velocity in terms 
of the osculating orbital elements: 



x = r cos cos (oj + f) — r cos i sin f2 sin (u + /), 
y = rsinSl cos (ui + f) + r cosi cos f2 sin (oj + /), 
z = rsinisin(w + /), 



(A5) 



and 



cos fi[r cos (ui + f) — rf sin (uj + /)] 

— cosisinil[fsin {uo + f) + rf cos (w + /)], 
sin fl[r cos (u> + /) — rf sin (w + /)] 

+ cos i cos f2[r sin (u + f ) + rf cos (uj + /)] , 
sinz[f sin (cj + /) + rf cos (w + /)], 

where r, f, and rf are in terms of a, e, and / (e.g., Murray & Dermott 1999). We also need 

dr r 



V = 



z = 



(A6) 



da 

dr 
~de 

df 
da 

dr 
de 

d(rf) 
da 

d(rf) 
de 



2er r 2 cos / 



a(l 



r 

Ya 



e(l-e 2 )' 



(A7) 



_rj_ 
2a 

r/(e + cos/) 
(1 -e 2 )(l + ecos/)' 



From equations (A3), (A5), and (A7), we find that 



dx 
~dt 



dx 
+ ~dt 



x . 

-a + 

a 



1 + e 2 



a(l-e 2 ) 1-e 2 



x . 
—e 

e 



(A8) 



the additional terms for dy/dt and dz/dt are similar. The additional terms for each of dx/dt, dy/dt, dz/dt 
are distinct for variations in e, and we have 



dx 
~dl 



dx 
~dl 



-—a + cosfi 
2a 



gcos(u, + /)-^ S in(u, + /) 
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dy 
dt 



dt 



dz 



dz 
~dt 



cos i sin £1 



sin 



y_ 

2a 

+ cos i cos f2 



2a 



-a + sin i 



^sin( W + /) + ^cos(u; + /) 



^cos(u; + /)-^sin(u; + /) 



gsin( W + /) + ^cos(u; + /) 



gsin(u, + /) + ^cos(u, + /) 



(A9) 



where equation (A7) should be used to get the functional forms. Unfortunately, the orbital elements must 
be calculated at each call to the differential equations when there is eccentricity damping. 
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Table 1 

Best Fit Orbital Parameters for the G J 876 Planets 



2-Kepler Fit a Dynamical Fit a 

Keck+Lick Keck Kcck+Lick 









sini = 


0.55 


sini = 


0.78 


Parameter 13 


Inner 


Outer 


Inner 


Outer 


Inner 


Outer 


M (Mj) 


0.56/ sini 


1.89/ sin i 


1.06 


3.39 


0.766 


2.403 


P (day) 


30.12 


61.02 


29.995 


62.092 


30.569 


60.128 


a (AU) 


0.130 


0.208 


0.1294 


0.2108 


0.1309 


0.2061 


e 


0.27 


0.10 


0.314 


0.051 


0.244 


0.039 


w (deg) 


330 


333 


51.8 


40.0 


159.1 


163.3 


T (JD) 


2450091.6 


2450602.09 


2449679.63 


M (deg) 


0.0 


-85.9 


289 


340 


356 


173 



a Two-Kcplcr fit by Marcy et al. 2001 and dynamical fit using a 
Levenberg-Marquardt iV-body integration scheme by Laughlin & Cham- 
bers 2001. 

b Thc parameters are the planetary mass M in terms of Jupiter mass 
Mj, the period P, the semi- major axis a, the orbital eccentricity e, the 
longitude of the periapse w, the epoch T, and the mean anomaly M. The 
stellar mass is 0.32M Q . 
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Table 2 

Comparison of contribution of various terms to the periapse precession rates for the 

Io-Europa AND GJ 876 SYSTEMS 



Io-Europa GJ 876 Planets 

Oi = 0°, 2 = 3 = 180° 0i = 2 = 03 = 0° 

d = 0.0026, e 2 = 0.0013 a = 0.255, e 2 = 0.035 



dzui/dt dwnjdt dw\jdt dzD2/dt 

Terms Ordcr a ^day" 1 ) ("day" 1 ) ("day" 1 ) ^day" 1 ) 



cos k9 3 


e 2 


0.0034 


0.0092 


0.0374 


-0.0470 




e 4 


0.0000 


0.0000 


0.0032 


-0.0080 


cos (fc0i - £0 2 ) 


e 1 


-1.4832 


-1.5778 


-0.2503 


0.1677 




e 2 


0.0190 


0.0820 


0.1456 


-0.3983 




e 3 


-0.0002 


-0.0011 


-0.0809 


0.2471 


Total 




-1.46 


-1.49 


-0.145 


-0.039 


Actual 










0.116 



a Terms of order e ^ +2m e)^ +2n ', with \k\ + \£\ + 2m + 2n = N, in the disturbing 
potential (eq. [8]) are grouped together under order e N . 
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